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Cundall& Stack (1979D or |Thornton (2000| ), except 
tha t stresses, rather than strai ns, are controlled, as 
in (|Parrinello & Rahman 198 Contact elasticity re- 
lates the normal force Fn is to the normal deflection 
h of the contact by the Hertz law 
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/i 3/2 , with E 
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while variations of tangential elastic forces Ft with 
tangential displacements 5ut, 
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O^BSTRACT: Elastic properties and internal states of isotropic sphere packings are studied by numerical sim- 
^lilations. Several numerical protocols to assemble dense configurations are compared. One, which imitates 
Experiments with lubricated contacts, produces well coordinated states, while another, mimicking the effect of 
i-^ibrations, results, for the same density, in a much smaller coordination number z, as small as in much looser 
ly-systems. Upon varying the confining pressure P, simulations show a very nearly reversible variation of density, 
^vhile z is irreversibly changed in a pressure cycle. Elastic moduli are shown to be mainly related to the coordi- 
nation number. Their P dependence notably departs from predictions of simple homogenization approaches in 
* ^he case of the shear moduli of poorly coordinated systems. 

Th ■ 

pi INTRODUCTION 

^jThe mechanical properties of solidlike granular mate- 
cjials are well known to depend on the internal struc- 
ture of the packing. Classically, one distinguishes 
—Ibetween the behaviour of dense and loose materi- 
als (I Wood 1990L Mitchell 19931). Howe ver, some re- 
Qults - see e.g. (Benahmed et al. 2004) - also indi- 

1 ^ate that other factors than the sole packing frac- 

tion (or void index) also determine the quasistatic re- 
sponse to applied load variations. In numerical sim- 
f-Oilations, it is a common practice to remove fric- 
OQion in the asse mbling stage in order to prepare 
£jense samples (IMakse et al. 19991 IThornton 20001) . 
Mft> is not guaranteed that the correct initial state, 
^^s obtained in laboratory experiments, is repro- 
Q<duced. Elastic properties, or sound wave veloci- 
Qies, are now commonly measured in soil mechan- 
ics (|Chen et al. 19881 IHicher 1996 ) and condensed 
^Tnatter physics dJia et al. 19991) laboratories. Their 
Vfvaluation in numerical calculations can allow for 

comparisons with experiments. 

We report here on a numerical study of isotrop- 

ically assembled and compressed sphere packings, 

prepared in different initial states. Coordination num- 
bers are found to vary according to the preparation 

method independently from density, and to determine 

the elastic moduli and their pressure dependence. 

2 MODEL AND NUMERICAL METHODS 

Numerical samples of 4000 identical balls of di- 
ameter a and mass m are prepared by stan- 
dard molecular dynamics calculations, involving 
periodic boundary conditions in all three di- 
rections. The method is similar to those of 



d¥ T 
d5u T 



dF, 



QLT' 



N 



dh 



with ot T 



2(1 



(2) 



are modeled with a simplified form of the Cattaneo- 
Mindlin-Deresiewicz theory (1 Johnson 198 5 ). On im- 
ple menting Q, special care was taken, as advocated 
by Elata & Berryman (1996), to avoid spurious cre- 
ation of elastic energy. 

Particles are endowed with the Young modulus 
E = 70 GPa and the Poisson coefficient v = 0.3 of 
glass beads. The tangential reaction is limited by the 
Coulomb condition with friction coefficient /x = 0.3. 

First, the sample is assembled from an initial 
disordered loose granular gas state, under the pre- 
scribed isotropic pressure P = 10 kPa. In proce- 
dure A, tangential forces are suppressed in this stage, 
as for frictionless grains. This produces dense sam- 
ples with a coordination number z*, counting only 
force-carrying grains and contacts, approaching 6 
in the rigid limit (|Roux 2000|) . Only a small frac- 
tion /o — 1.3% of grains (the "rattlers") carry no 
force. This procedure, already used in other numer- 
ical work (IThornton 20001) . amounts to dealing with 
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perfectly lubricated beads. Imperfect lubrication can 
be modelled with a very small friction coefficient, 
/io = 0.02, resulting in slightly different samples, de- 
noted as B configurations. 

A different assembling procedure (method C) was 
designed to simulate dense samples obtained by vi- 
bration. Configurations A are first dilated, scaling all 
coordinates by a common factor A = 1.005, thereby 
supressing the contacts ; then the grains are attributed 
random velocities and mixed with a kinetic-energy 
preserving event-driven ("hard sphere") calculation, 
until each of them has undergone 50 collisions on av- 
erage ; finally, they are compressed with strongly dis- 
sipative collisions and friction, to mechanical equilib- 
rium under P = 10 kPa. Remarkably (see table [I]), 
such configurations C are very nearly as dense as the 
perfectly lubricated ones A, and actually denser than 
B ones, while their active coordination number z* is 
considerably lower, with many more rattlers. 

Finally, procedure D consists in directly compress- 
ing the granular gas to equilibrium at 10 kPa with the 
final coefficient of friction ji = 0.3. z* values are close 
to the C case, but the density is significantly lower. 

Tabled] summarizes the data on these initial states. 
All data throughout this paper are averaged over 5 
different samples, error bars correspond to one r.m.s. 
deviation. These initial states are then further corn- 
Table 1. Packing fraction <t>, coordination number z* on force- 
carrying structure and proportion of rattlers fo at the lowest pres- 
sure 10 kPa for the four simulated preparation procedures. 



exhibit very similar properties at high P. The shape 



State 




z* 


fo (%) 


A 


0.637 ±0.009 


6.074 ±0.002 


1.3 ±0.2 


B 


0.627±2-10~ 4 


5.80 ±0.007 


1.65 ±0.02 


C 


0.635 ±0.002 


4.56 ±0.03 


13.3 ±0.5 


D 


0.606 ±0.002 


4.62 ±0.01 


10.4 ±0.9 



pressed, on applying pressure steps, up to 100 MPa, 
assuming contacts still behave elastically, and then the 
pressure is gradually decreased back to 10 kPa. To en- 
sure quasistatic conditions are maintained with suffi- 
cient accuracy, strain rates e are constrained by con- 
dition t-\p^p < 10~ 4 . Equilibrium states are recorded 

for pressure values at ratio vTO. Throughout this pro- 
cess, the friction coefficient is maintained equal to 
0.3, for all four configuration types. Elastic constants 
are measured on building the stiffness matrix associ- 
ated with the contact network at equilibrium. 

3 RESULTS 

3 . 1 Structure of equilibrium configurations 
The results obtained on samples A, C and D are re- 
ported below, configurations of type B behaving very 
similarly to A ones. Fig. Q] displays the evolution of 
packing fraction $ and coordination number z* over 
the compression-decompression cycle. Changes in $ 
are very nearly reversible (elastic): density differences 
between states A, C and D, are maintained at low 
pressure after one cycle, although A and C samples 
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Figure 1. Variations of $ and z* as P increases up to 100 MPa 
and decreases back to 10 kPa, in samples A (square dots), C 
(open circles), and D (crosses). 



of the normal force distribution also changes with P. 
It can be characterized by the reduced moments: 
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while the average normal force, for monosized beads, 
simply relates to P as 
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z = z*(l — fo) being the total coordination number. 
The width of this distribution, as expressed, e.g., by 
Z(2), decreases at growing P, the fastest in well- 
coordinated A samples. Changes in friction mobiliza- 
tion are also observed : as P grows, it first decreases 
in C and D samples, and increases in A ones (in which 
it starts at zero). 

The change of z* in type A samples as P decreases 
from a high value - many more contacts are lost than 
were gained at increasing P - might seem surprising. 
One should note however that configurations with a 
high coordination number, for nearly rigid grains, are 
extremely rare. Each contact requires a new equation 
to be satisfied by the set of sphere centre positions. 
Equilibrium states of rigid, frictionless sphere assem- 
blies, apart from the motion of the scarce rattlers, are 
isolated points in configuration space, because of iso- 
staticity (Roux 2000). As isotropic compression, at 
the microscopic scale, is not reversible, due to friction 
and to geometric changes, one should not expect ex- 
ceptional configurations to be retrieved upon decreas- 
ing the pressure. Large coordination numbers of A (or 
B) samples do not survive a pressure cycle. The his- 
tory of an isotropic sample can therefore significantly 
influence its structure without any appreciable density 
change. 
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3.2 Elastic moduli 

Bulk (B) and shear (G) moduli of all equilibrium 
states for ascending P were computed on solving 
linear systems of equations involving the stiffness 
matrix (they express the response to infinitesimal 
stress changes). Their variations with P are plotted 
on Fig. [2] Results for configurations B (not shown) 
are very close to those of A samples. The obvious first 
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Figure 2. P dependence of bulk moduli B (bottom) and shear 
moduli G (top) (symbols for A, C, D states as on fig. [T] joined 
by continuous lines), and of their upper (B and G) and lower (for 
B only) bounds for states A and C (same symbols, no line). Note 
the relatively narrow bracketing of B in all cases, and the large 
overestimation of G by its Voigt upper bound in C samples. 

conclusion to be drawn is that elastic moduli are sen- 
sitive to coordination rather than density, as results for 
states C and D are very similar. The pressure depen- 
dence of bulk moduli differs a little between A sam- 
ples on the one hand and C, D on the other. In the 
latter case, the increases of B with P is slightly faster 
than the P 1 / 3 law predicted by simple estimates (see 
below). The most striking behaviour is that of G in 
samples C and D, the increase of which approaches a 
P 1 / 2 dependence. To explain such observations, one 
can try to estimate the moduli as follows. Assuming 
the distribution of normal forces is known, one gets 
by virtue of (OQ) and © the distribution of contact 
stiffnesses. It is easy, then, to derive upper bounds 
to B and G, and a lower bound to B, analogous to 
the Voigt and Reuss bounds for elastic heterogeneous 
continua (Nemat -Nasser & Hori 19931 . The Voigt up- 
per bound is the simple "effective medium" estimate 
that results from the assumption of affine displace- 



ment fields. Using the properties and notations intro- 
duced in Eqns.CD [HE and |H one gets: 
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(5) 



-Q Voigt 



To write a lower bound (Reuss estimate), one needs a 
trial set of equilibrium contact force increments cor- 
responding to the stress increment. For a simple in- 
crease of isotropic pressure, this is readily obtained 
by a scaling of the forces corresponding to the pre- 
existing pressure. Hence a lower bound for B (but no 
such estimate is available for G). Denoting as ttn the 
ratio 1 1 Ft | |/P/v in each contact, and defining 
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a modified reduced moment Z(5/3) (Eqn.O, one has: 



B > B 



Reuss 




(6) 



In view of the force distribution and mobilization of 
friction observed, B, bracketed by © and [6l cannot 
depart very much from a P 1 / 3 dependence {r\ N jctT 
is at most 0.11 anyway f or [i — v = 0.3, and the prod- 
uct Z(l/3)Z(5/3) only exceeds 1.15 for systems with 
few contacts at low pressure). As to the increase of z 
with P, it does not appear to entail large effects either. 

The behaviour of G for samples C and D is quite 
different. G seems to get unexpectedly small (G < 
B/2) at low pressure. Its upper bound is a very 
poor estimate in such cases (as, by ©, G Ymgt = 
1.34 x P Voigt ). The situation is reminiscent of friction- 
less sphere packings (O'He rnet al. 20031) . for which 
G <^ B under isotropic pressure in the rigid limit. We 
could observe that those states had the largest level 
of strain fluctuations (departure from affine displace- 
ment field). 

4 CONCLUSION 

Our simulations of isotropically assembled sphere 
packings revealed the following points. 

• Configurations of a given density can vary con- 
siderably in coordination number. Samples as- 
sembled with a procedure designed to imitate vi- 
bration can have a large density and a small co- 
ordination number. 

• Elastic moduli are primarily sensitive to coordi- 
nation numbers. 
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• They vary with pressure in rather good agree- 
ment to simple predictions (nearly as P 1 / 3 , with 
a small effect of contact creation as P increases) 
in highly coordinated samples, but the shear 
modulus behaves quite anomalously in poorly 
coordinated ones, with a low value at low P and 
a faster increase, nearly as P 1 / 2 in some cases. 

• A compression-decompression cycle, although 
almost reversible in terms of density, can sub- 
stantially reduce the coordination number when 
it was initially high. 

We therefore suggest to use elastic moduli, which can 
be compared between simulations and experiments, 
as indicators of the internal state (contact density) of 
granular packings. 

On comparing numerically predicted ultrasonic 
wave velocities with experimental values obtained on 
dense sphere packs with P in the range 100-800 kPa, 



we observe ( Agnolin et al. 2005 1 that perfectly lubri 



cated samples (type A or B), are considerably too 
stiff, even though they agree with experimental ob- 
servations in the MPa range (Maks e et al. 19 99). Al- 
though somewhat idealized, our "vibrated" ones (type 
C) seem to be closer to the materials studied in the 
laboratory. We also note in another contribution to 
the present proceedings (IRoux 20 05 ) that their stress- 
strain curves under growing deviatoric stress are also 
closer to experimentally observed mechanical be- 
haviours. Of course, it will be necessary in the near fu- 
ture to investigate by numerical simulations more "re- 
alistic" assembling procedures (Ema m et al. 20 05 ). 
and the effects of the resulting anisotropy. 
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